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ABSTRACT 

We present a study of the extent to which the Sehwood & Binney radial migration of 
stars is affected by their vertical motion about the midplane. We use both controlled 
simulations in which only a single spiral mode is excited, as well as slightly more 
realistic cases with multiple spiral patterns and a bar. We find that rms angular 
momentum changes are reduced by vertical motion, but rather gradually, and the 
maximum changes are almost as large for thick disc stars as for those in a thin disc. 
We find that particles in simulations in which a bar forms suffer slightly larger angular 
momentum changes than in comparable cases with no bar, but the cumulative effect 
of multiple spiral events still dominates. We have determined that vertical action, and 
not vertical energy, is conserved on average during radial migration. 

Key words: galaxies: kinematics and dynamics - galaxies: evolution - galaxies: 
structure. 



1 INTRODUCTION 



It has recently become clear that galaxy discs evolve signifi- 
cantly over time due to both internal and external influences 
(see van der Kruit Freeman]|2011 Sellwood]|2010 for re- 



views). Scattering by giant molecular clouds {e.g. Spitzer 
|& Schwarzschild 1953) is now thought to play a minor role 



([Hanninen & Flynn 2002J, while spirals (Sellwood & Carl- 


berg|1984 


Binney & Lacey||1988 Sellwood & Binney 


2002) 


and interactions (Helmi et a/.|2006 Kazantzidis et al. 


2008) 



are believed to cause more substantial changes. 

The Milky Way is believed to have both a thin and a 
thick disc ([Gilmore &: Reid|1983| [Munn et a/.|2004| |Juric~ 



a/.|2008| Ivezic et a/.|2008" ), as do perhaps many disc galax- 
ies ( Yoachim Dalcanton|[200 6). Relative to the thin disc, 
the thick disc has a higher velocity dispersion, lags in its 
net rotational velocity ( Chiba Beers|2 000), contains older 
stars with lower metallicities (Majewski 1993|) and enhanced 
[a/Fe] ratios (Bensby et al. 2005; Red dy, Lambert &: Allende 
Prieto|2006l |Fuhrmann|2Q08| . A number of different criteria 
have been used to divide stars into the two disc populations, 
and some of these trends may depend somewhat on whether 
a spatial, kinematic, or chemical abundance criterion is ap- 



plied (Fuhrmann 2008 Schonrich & Binney 2009b Lee et 



al. 2011). Recently, Bovy et al. (2011) suggest that there 



is no sharp distinction into two populations, but rather a 
continuous variation in these properties. 

Several models have been proposed for the formation 
of the thick disc: accretion of disrupted satellite galaxies 



(Abadi et al. 2003), thickening of an early thin disc by a 



|bos Helmi||2008 ), star formation triggered during galaxy 



assembly (Brook et a/.[2005 Bournaud, Elmegreen & Mar- 



itig,2009t, and radial migration ( Schonrich Binney|2009a 
Loebman et a/.|2011 ) . While it is likely that no single mech- 



anism plays a unique role in its formation, our focus in this 
paper is on radial migration. 



minor merger (Quinn, Hernquist & Fullagar 1993| Villalo- 



As first shown by Sellwood & Binney ( 2002 ) , spiral ac- 
tivity in a disc causes stars to diffuse radially over time. 
Stars near corotation of a spiral pattern experience large 
angular momentum changes of either sign that move them 
to new radii without adding random motion. As spirals are 
recurring transient disturbances that have corotation radii 
scattered over a wide swath of the disc, the net effect is that 
stars execute a random walk in radius with a step size rang- 
ing up to 2 kpc, which results in strong radial mixing. 
Note that this mechanism is distinct from the "blurring" 
caused by epicycle oscillations of stars about their home 
radii, Rh, where LI = Rl\d^/dR\ in the midplane and 
$ is the gravitational potential. Radial mixing is dubbed 
"churning" and causes the home radii themselves to change. 

Note that radial mixing caused by recurrent transient 
spirals has two related properties that distinguish it from 
other processes that have sometimes also been said to cause 
radial migration. The two distinctive properties are the ab- 
sence of disc heating, already noted, and the fact that stars 
mostly change places, causing little change to the distribu- 
tion of angular momentum within the disc and leaving the 
surface density profile unchanged. These are both features 
of scattering at corotation; spirals also cause smaller angular 
momentum changes at Lindblad resonances which do heat 
and spread the disc. 
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Bar formation causes greater angular momentum 



changes than those of an individual spiral pattern (Friedli, 
Benz &; Kennicu tt 1994; Rab oud a/.||1998| [Gr enon 1999; 



Debattista et a/.i2QQ6; MinchevTt a/.|2011||Bird, Kazantzidis 



&: Weinberg||2011| [Brunetti et a/.||2Qll ). However, because 
the bar persists after its formation, the associated angular 
momentum changes have a different character from those 
caused by a spiral that grows and decays quickly. Bar for- 
mation changes the radial distribution of angular momen- 
tum, and therefore also the surface density profile of the 
disc, as well as adding a substantial amount of radial mo- 
tion ( |Hohl||1971| . Although the bar may subsequently set- 
tle, grow, and/or slow down over time, the main changes 
associated with its formation probably occur only once in 
a galaxy's lifetime (for a dissenting view see Bournaud 
Combes||2002 ) whereas spirals seem to recur indefinitely as 



long as gas and star formation sustain the responsiveness 
of the disc. Thus angular momentum changes in the outer 
disc beyond the bar should ultimately be dominated by 
the effects of radial mixing by dynamically-unrelated tran- 
sient spirals, as we find in this paper. Note that Minchev fc| 
[Fam aey (2010) studied the separate physical process of reso- 
nance overlap using steadily rotating potentials to represent 
a bar and spiral, and therefore did not address possible mix- 
ing by transient spirals in a barred galaxy model. 

f2009) and Bird, Kazantzidis & 



Both Quillen et al. 



[Weinberg ( 2011 ) find that external bombardment of the disc 



can enhance radial motion through the increase in the cen- 
tral attraction and also through possible angular momen- 
tum changes. However, the importance of bombardment is 
unclear since satellites may be dissolved by tidal shocking 
( D'Onghia et a/.||2010 |) before they settle, and furthermore 
satellite accretion to the inner Milky Way over the past 
~ 10^° years is strongly constrained by the dominant old 
age of thick disc stars {e.g. [Wyse 2009). 



Sellwood & Binney (2002| discovered spiral churning in 



2D simulations, and it has been shown to occur in fully 3D 
simulations (Roskar et al. 2008a|b Loebman et al. 2011). 
Schonrich & Binney (2009a) present a model for Galactic 



chemical evolution that includes radial mixing, and point 



out ( Schonrich Binney||2009a|b Scannapieco et a/.||2011 ) 
that it naturally gives rise to both a thin and a thick disc, un- 
der the assumption that thick disc stars experience a similar 
radial churning. Loebman et al. (2011) show that extensive 



spiral activity in their simulations causes a thick disc to de- 



velop, and present a detailed comparison with data (Ivezic et 
ari|2008t from SDSS ([York et a/.||2000|. "Sanchez-Blazquez 



et al. (2009) and [Martmez-Serrano et al. (2009| also gen- 



erated galaxies with realistic break radii in the exponen- 
tial surface brightness profiles using cosmological A/'-body 
hydrodynamic simulations that included radial migration. 



However, Bird, Kazantzidis & Weinberg (2011) report that 



no significant thick discs developed in their simulations. 



Lee et al. (2011) find evidence for radial mixing in the 



thin disc and for a downtrend of rotational velocity with 
metallicity. The presence of such a trend does not rule out 
migration in the thick disc, but shows that mixing for the 
oldest stars cannot be complete. Bovy et al. (2011) find a 



continuous change of abundance-dependent disc structure 
with increasing scale height and decreasing scale length 
which they note is the "almost inevitable consequence of 
radial migration" . The decreasing metallicity gradient with 



age of Milky Way disc stars can also be attributed to radial 

mixing ( |Yu et ai.||2012t. 

While Loebman et al. (2011 ) present circumstantial evi- 



dence for radial migration in the thick disc, they do not show 
explicitly that it is occurring, neither do they attempt to 
quantify the extent to which it may be reduced by the weak- 
ened responses of thick disc stars to spiral waves in the thin 
disc. Bird, Kazantzidis & Weinberg" f2011") find that mixing 
is more extensive when spiral activity is invigorated by star 
formation, although the level of spiral activity is strongly de- 
pendent on the "gastrophysical" prescription adopted. They 
show that mixing persists even for particles with large oscil- 
lations about the mid-plane, and they determine migration 
probabilities from their simulations. 

In this paper, we set ourselves the limited goal of deter- 
mining the extent of radial migration in isolated, collisionless 
discs with various thicknesses and radial velocity dispersions 
that are subject to transient spiral perturbations. Following 



Sellwood & Binney (2002), we first present controlled simu- 



lations of two-component discs constructed so as to support 
an isolated, large-amplitude spiral wave, in order to study 
the detailed mechanism of migration in the separate thin and 
thick discs. We also report the responses of two-component 
discs to multiple spiral patterns, both with and without a 
bar. 



2 DESCRIPTION OF THE SIMULATIONS 

All our models use the constant velocity disc known as the 
Mestel disc f Binne y Tremaine|2 008 , hereafter BT08). Lin- 
ear stability analysis by Toomre ( ^1981^ see also Zang 1976 
and Evans & Read 1998) revealed that this disc with moder- 
ate random motion lacks any global instabilities whatsoever 
when half its mass is held rigid and the centre of the remain- 
der is cut out with a sufficiently gentle taper We therefore 
adopt this stable model for the thin disc, and superimpose 
a thick disc of active particles that has 10% of the mass of 
the thin disc. The remaining mass is in the form of a rigid 
halo, set up to ensure that the total central attraction in 
the midplane is that of the razor-thin, full-mass, untapered 
Mestel disc. 



2.1 Disc set up 

Ideally, we would like to select particles from a distribution 
function (DF) for a thickened Mestel disc. Toomre (1982) 
found a family of flattened models that are generalizations 
of the razor-thin Zang discs that have the isothermal sech^ 



vertical density profile ( |Spitzer|[l942 Cammp950 ). Unfor- 
tunately, these two-integral models have equal velocity dis- 
persions in the radial direction and normal to the disc plane 
(see BT08), whereas we would like to set up models with 
flattened velocity ellipsoids. 

Since no three-integral DF for a realistic disc galaxy 
model is known, as far as we are aware, we start from the 
two-integral DF for the razor- thin disc ( |Zang||1976| |Toomre| 
T977k: 



^ I Sellwood] ( |2012| finds that particle realizations of this disc are 
not completely stable, but on a long time scale when N is large. 
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(1) 



where ^ is a particle's specific energy and Lz its specific z- 
angular momentum. The free parameter q is related to the 
nominal radial component of the velocity dispersion through 

an = Voil + qr''\ 



(2) 



with Vb being the circular orbital speed at all radii. The 
value of Toomre's local stability parameter for a single com- 
ponent, razor-thin Mestel disc would be 



Q 



3.36/(1 + 5)1/2' 



(3) 



where / is the active fraction of mass in the component; note 
that these expressions for both gr and Q are independent 
of radius. A composite model having two thickened discs, 
such as we employ here, will have some effective Q that is 
not so easily expressed. We choose different values of gr for 
the thin and thick discs, given in Table [l] below, in order 
that the thick disc has a greater radial velocity dispersion, 
as observed for the Milky Way. 

We limit the radial extent of the discs by inner and 
outer tapers 



ME,L.) = 



[1 + (L,/L,)4][1 + (L,/L„)6]' 



(4) 



where Li and Lo are the central angular momentum values of 
the inner and outer tapers respectively. The exponents are 
chosen so as not to provoke instabilities (Toomre, private 
communication). We employ the same taper function for 
both discs and choose Lo = 15L^. We further restrict the 
extent of the disc by eliminating all particles whose orbits 
would take them beyond 25Ri, where Ri = Li/Vo is the 
central radius of the inner taper. This additional truncation 
is sufficiently far out that the active mass density is already 
substantially reduced by the outer taper. 

We thicken these disc models by giving each the scaled 
vertical density profile 

p(z) oc ^, (5) 

^ ^ (el-~l/2+0.2e-5|-~l/2)2' ^ ^ 

where z — z/zq with zq independent of R. While both this 
function and the usual sech^z function have harmonic cores, 
we prefer eq. ([s]) because it approaches p oc e~'^' more 



rapidly, as suggested by data {e.g. Kregel, van der Kruit & 
de Grijs||2002 ); zq is therefore the exponential scale height 
when z ^ zq. We set the vertical scale height of the thick 
disc to be three times that of the thin disc, as suggested for 
the Milky Way ( (Juric et a/.||2008| . 

We estimate the equilibrium vertical velocity disper- 
sion at each z-height by integrating the ID Jeans equation 
(BT08, eq. 4.271) in our numerically-determined potential. 
This procedure is adequate when the radial velocity disper- 
sion is low, but the vertical balance degrades in populations 
having larger initial radial motions, which fiare outwards as 
the model relaxes from initial conditions, as shown below. 

The tapered thin Mestel disc has a surface density 
that declines with radius as shown in Fig. ^ While there 
is no radial range that is closely exponential, we estimate 
an approximately equivalent exponential radial scale length 
Rd — 4.0i?i as shown by the blue line. Accordingly, we 
choose zq = O.ARi for the thin disc so that the ratio Rd/zQ 
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Figure 1. Initial surface density profiles of the tapered thin (top 
black curves) and thick (bottom green curves) Mestel discs in sim- 
ulation M2 (solid curves) measured from the particles. The dashed 
red curves, which are almost perfectly overlaid by the solid curves, 
show the expected surface density of the tapered discs from inte- 
grating the DP over all velocities, while the dotted curves indicate 
the density profiles of the corresponding non-tapered discs. The 
groove in the thin disc is centred on = Q.hRi and is broadened 
by epicycle excursions. The blue line of slope —0.25 indicates the 
adopted exponential profile of the thin disc. 



is similar to that of the Milky Way and the average volume- 
corrected ratio of 7.3 =b 2.2(lcr) found by [Kregel, van der 
Kruit & de Grijs ( 2002 ) for 34 nearby galaxies. The values 



we adopt for each model are given in Table [l] below. 

The radial gravitational potential gradient in the mid- 
plane of our model is less than that of the full-mass, razor- 
thin, infinite Mestel disc as a result of the reduced disc sur- 
face density, the angular momentum tapers, the finite thick- 
ness of the discs, and gravitational force softening. In order 
to create an approximate equilibrium, we therefore add a 
rigid central attraction to the self- consistent forces from the 
particles in the discs at each step. We tabulate the supple- 
mentary central attraction needed at the initial moment to 
yield a radial force per unit mass of —Vq/R in the midplane, 
and apply this unchanging extra term as a spherically sym- 
metric, rigid central attraction. 

We find that our model is close to equilibrium, with an 
initial virial ratio of the particles of ~ 0.52. This value ad- 
justs quickly to reach a steady virial ratio of 0.50 within the 
first 32 dynamical times (defined below). As noted above, 
the initial imbalance seems to arise mostly from the verti- 
cal velocity structure, for which we adopted the ID Jeans 
equation. The adjustment of the model to equilibrium is 
illustrated in Fig. [2] the bottom left panel shows that thick- 
nesses of both discs change as the system relaxes, increasing 
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Table 1. Parameters of all the simulations described in this paper. 



Simulation 


Disc 




/ 


Ri 


Vo 


N 


m 


Az 


£ 

Ri 


Rings X Spokes x Planes 


M2 


thin 




0.5 


0.4 


0.283 


1200 000 


, 2 


0.1 


0.15 


120x128x243 




thick 




0.05 


1.2 


0.567 


1 800 000 










T 


thin 




0.5 


0.4 


0.283 


1200 000 





0.1 


0.15 


120x128x243 




thick 




0.05 


1.2 


0.567 


1800 000 










M3 


thin 




0.3333 


0.4 


0.189 


1 200 000 


, 3 


0.1 


0.15 


120x128x243 




thick 




0.0333 


1.2 


0.378 


1800 000 










M4 


thin 




0.25 


0.4 


0.142 


1200 000 


, 4 


0.1 


0.15 


120x128x243 




thick 




0.025 


1.2 


0.283 


1800 000 










M4b 


thin 




0.25 


0.2 


0.142 


1200 000 


, 4 


0.05 


0.075 


120x128x405 




thick 




0.025 


0.6 


0.283 


1 800 000 












massless thick 





1.2 


0.283 


1800 000 










M2b 


thin 




0.5 


0.4 


0.283 


480 000 


, 2 


0.1 


0.15 


75x80x625 




thick 




0.05 


1.2 


0.567 


240 000 












massless 


1 





0.5 


0.567 


240 000 












massless 


2 





0.6 


0.567 


240 000 












massless 


3 





0.7 


0.567 


240 000 












massless 


4 





0.8 


0.567 


240 000 












massless 


5 





1.6 


0.567 


240 000 












massless 


6 





2.0 


0.567 


240 000 












massless 


7 





2.4 


0.567 


240 000 










M2c 


thin 




0.5 


0.4 


0.283 


480 000 


, 2 


0.1 


0.15 


75x80x243 




thick 




0.05 


1.2 


0.567 


240 000 












massless 


1 





1.2 


0.378 


240 000 












massless 


2 





1.2 


0.454 


240 000 












massless 


3 





1.2 


0.680 


240 000 












massless 


4 





1.2 


0.794 


240 000 












massless 


5 





1.2 


0.907 


240 000 












massless 


6 





1.2 


1.021 


240 000 












massless 


7 





1.2 


1.134 


240 000 










TK 


thick 




0.5 


1.2 


0.283 


1800 000 


, 2 


0.1 


0.15 


120x128x243 


UC , UCBl, 


thin 




0.4 


0.4 


0.227 


200 000 


^ 8, / 1 


0.2 


0.3 


75x80x125 


UCB2 


thick 




0.04 


1.2 


0.454 


300 000 











The first column gives the simulation designation used in the text. The next five columns give the 
properties of each disc, one for each line: / is the mass fraction, zq its vertical scale height, aji 
is the nominal radial velocity dispersion, and is the number of particles in the disc. The final 
four columns give m, the active sectoral harmonic(s), Az the vertical spacing of the grid planes, 
£ the softening length, and the grid size. The "massless" discs of simulations M4b, M2b, and M2c 
describe massless thick discs composed of test particles, zq, Az, and e are in terms of Ri, and aji 
is in terms of the circular velocity Ve- 



in the outer parts and decreasing in the inner parts - notice 
that the change in thickness is least over the radius range 
2Ri ^ R ^ lORi, where the surface mass density is closest to 
its untapered value (dotted curves in Fig.[T]). The changes, 
which are larger for the radially-hotter thick disc, result from 
the radial excursions of the particles, which may take them 
far from their initial radii for which the vertical velocity 
was set. However, neither the radial balance (top left panel) 
nor the vertical velocity dispersion (right panel) of the thick 
disc change significantly from their initial values. After this 
initial relaxation from initial conditions, we do not observe 
any significant flaring. A negligibly small fraction (< 1%) of 
particles escape from our grid by the end of the simulation, 
and most particle loss takes place as the model settles. 



Following Sellwood & Binney (2002), we seed a vigor- 
ously unstable spiral mode by adding a Lorentzian groove 
in angular momentum to the DP of the thin disc only: 



f{E,L,) = ME, L. 



1 + 



(6) 



Here /3, a negative quantity, is the depth of the groove, 
iL^L is its width, and is the angular momentum of the 
groove centre. This change to the DP seeds a predictable 
spiral instability ( [Sellwood &: Kahn||1991[ ). The groove in 
the thin disc has the parameters (3 — —0.9, wl — 0.3i?iVb, 
and L* = 6.5i?iVb. We find that a deeper and wider groove 



is needed than that used by Sellwood & Binney (2002) in 
order to excite a strong spiral in our thickened disc. 

Since a groove of this kind will provoke instabilities at 
many sectoral harmonics, m, we restrict disturbance forces 
from the particles to a single non-axisymmetric sectoral har- 
monic to prevent other modes from growing. Corotation for a 
global spiral mode excited by this groove is at a radius some- 
what greater than L*/Vb, where local theory would predict. 
A large-scale mode is not symmetric about the groove centre 
due to the geometric variation of surface area with radius 
- an effect that is neglected in local theory. However, the 
radius of corotation approaches the local theory prediction 
for modes of smaller spatial scale, z.e. as m ^ oo. 

As in iSellwood & Binney ( 2002 ) , we choose the circu- 
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Figure 2. The circular velocity in disc midplane (top left) and 
vertical velocity dispersion of the thick disc particles (right) at five 
equally spaced times during simulation M2. The times, which in- 
clude the initial and final moments of the simulation, are colour 
coded in temporal order by solid black, red, green, blue, and dot- 
ted cyan respectively. The dashed horizontal line in the top left 
panel shows the theoretical circular velocity Vq = 1. The bottom 

left panel shows the initial ^2;^^ of the thin (bottom curves) 
and thick (top curves) discs (solid black) compared to that at 
t = 32Ri/Vo (dashed red). 



lar velocity Vb to be our unit of velocity, Ri our unit of 
length, the time unit or dynamical time is tq = Ri/Vo, our 
mass unit is Mo = VqRi/G^ and Li is our unit of angu- 
lar momentum. One possible scaling to physical units is to 
choose Ri — 0.75 kpc, with Rd — 4:.0Ri being the equivalent 
scale length of the thin disc, and tq = 3.0 Myr, leading to 
Vb = 244 km s"^ and Mo = 1.04 x 10^° M©. 



2.2 Numerical procedure 



We use the 3D polar grid described in [Sellwood Val-] 
luri (1997) to determine the gravitational field of the par- 



ticles. This "old-fashioned" method is not only well suited 
to the problem, but also has the advantage of being tens of 
times faster than the "modern" methods recently reviewed 
by Dehnen & Read (2011). We employ a grid having 120 
rings, 128 spokes, and 243 vertical planes, and adopt the cu- 
bic spline softening rule recommended by Monaghan] ( |1992| ) , 
which yields the full attraction of a point mass at distances 
greater than two softening lengths (2^). The Plummer rule 
used in [Sellwood Binney ( 2002| is suited for razor-thin 
discs where it mimics the effect of disc thickness but, because 
it weakens inter-particle forces on all scales, it is unsuited 
for 3D simulations where disc thickness is already included 
in the particle distribution. The value of e, given in Table [l] 
exceeds the vertical grid spacing in order to minimize the 
grid dependence of the inter-particle forces. 

In the simulations described in §§3&4, we use quiet 



starts (Sellwood & Athanassoula 1986) for both discs to 



reduce the initial amplitude of the seeded unstable spiral 
mode far below that expected from shot noise. In 3D, this 
requires many image particles for each fundamental parti- 
cle: we place two at each (i?, 0, z) position with oppositely 
directed 2;- velocities and two more reflected about the mid- 
plane at the point (i?, 0, — z). We then place images of these 
four particles at equal intervals in 0, each set having identi- 
cal velocity components in cylindrical polar coordinates. For 



these simulations, the thin disc has 1200 000 particles and 
the thick disc 1800 000. All the particles in a single popu- 
lation have equal mass, but particle masses differ between 
populations in order to create the desired ratio of disc sur- 
face densities. 

We evaluate forces from the particles at intervals of 
0.02ro, and step forward some of the particles at each eval- 
uation. Since the orbital periods of particles span a wide 
range, we integrate their motion when R > 2Ri using longer 
time steps in a series of flve zones with the step doubling in 
length for each factor 2 in radius ( Sellwood||1985 ). We also 
subdivide the time step for particles within R = O.bRi with- 
out updating forces, with further decreases by a factor of 2 
for every factor of 2 decrease in radius (Shen & Sellwood 



2004). Note that very few particles have R < 0.5Ri, which 



is well within the inner taper where the rigid component of 
the force dominates. Tests with shorter time steps yielded 
similar results. 



2.3 List of simulations 

For convenience, we summarize the parameters of all the 
simulations presented in this paper in Table ^ Simulation 
T is constrained to remain axisymmetric, while all those 
whose identifler begins with 'M' are highly controlled exper- 
iments designed to support a single spiral instability. We first 
present, in §3, a detailed description of M2, which supports 
a bisymmetric spiral, and briefly compare it to simulation T. 
Variants of M2, with many populations of test particles are 
presented in §3.4, while simulation TK, which has a thick 
disc only, is described in §3.5. Simulations that support a 
single spiral of higher angular periodicity are motivated and 
described in §4. The last three simulations in the Table, with 
identifler beginning with 'U', are uncontrolled experiments 
presented in §5 that explore the consequences of multiple 
spirals, with and without a bar. 



3 A SINGLE BISYMMETRIC SPIRAL 

Our flrst objective is to study radial migration in both the 
thick and thin discs due to a single spiral disturbance. We 
therefore present simulation M2, which is designed to sup- 
port an isolated m — 2 spiral instability. 
We measure 
rR2 

A^{t) ^271 S(i?, 0, t)e^"^'^ dR, (7) 

Jri 

where ll{R^(j)^t) is the vertically-integrated mass surface 
density of the particles at time t, and we generally choose 
R^ = 1.5i?i and R2 = 19i?^. 

The top panel of Fig. [3] shows the time evolution of 
A2/A0, revealing that the mode grows exponentially, peaks, 
and then decays non-exponent ially to a trough. The peak is 
2.55 times higher than the minimum of the following trough. 
Continued evolution reveals that the amplitude rises again 
due to the growth of a secondary wave. In order to isolate the 
single initial spiral, we stop the simulation at the trough and 
compare quantities, such as the speciflc angular momenta of 
particles, at this flnal time with those at the initial time t — 
O.ORi/Vo. The bottom panel plots the same quantity from 
simulation T in which disturbance forces from the particles 
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Figure 3. The time evolution of A2/A0 in simulations M2 (top) 
and T (bottom). Note the difference in the vertical scales. The 
numbered arrows in the top plot mark the six times of M2 at 
which the snapshots of the discs are shown in Fig. [4] The top axis 
shows time scaled to physical units using the adopted scaling at 
the end of section 2.1. 
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were restricted to the axisymmetric {m — 0) term; the very 
slow rise is caused by the gradual degradation of the quiet 
start. 

Fig.[4]shows snapshots of the thin and thick discs of M2 
at six different times of the simulation that are marked by 
the red arrows in Fig. [3] 

We find the period of exponential growth can be very 
well fitted by a single mode, using the apparatus described in 
Sellwood & Athanassoula (1986). Table [2] reports the fitted 
pattern speed and growth rate 7, together with other 
parameters of the spiral in M2. Furthermore, the no n- linear 
evolution visible in Fig. [4] shows no evidence of a bar. Thus, 
the selected time period of the simulation M2 presents the 
opportunity to study radial migration due to a single, well 
isolated spiral wave. Note that the spiral pattern makes a 
full rotation every 27r/Qp :^ 45. Si^^/Vb. 



3.1 Angular momentum changes 

Fig. [5] shows the change in the specific 2;-angular momenta 
AL2 of the particles in the thin disc (top) and the thick 
disc (middle) of simulation M2 against their initial Lz. The 
deficiency of particles at Lz{t = 0) = = 6.5i?iVb in the 
top panel is due to the groove in the thin disc. The vertical 
lines show the locations of the corotation (solid blue) and 
the Lindblad resonances (dashed green) for nearly circular 
orbits. 

The maximum changes in angular momentum occur 
near corotation, and lie close to the solid blues line of 
slope —2; these particles cross corotation, in both direc- 
tions, to about the same radial distance away from it as 
they were initially. As for the razor-thin disc, we find that 
large angular momentum changes occur only around the 
time that the spiral saturates, for reasons explained by 
Sellwood fc^illlliyl ( |2002t . While many particles near the 





> 



(N 



> 




> 






Figure 4. Snapshots of the thin (left column) and thick (right 
column) discs in simulation M2 at the six times marked in Fig. [3] 
One in every 24 particles is plotted for the thin disc, and one in 
every 36 for the thick. The side views are shown only at the initial 
time because they change little throughout the simulation. 
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Table 2. Measured values from fits to the spiral modes. 





Tfl 


O / ^0 

^^P' Ri 


7 


Rc 
Ri 


(Am,) 

V Ao / peak 


P /T 


vq ' v^^y^j 


M2 


2 


0.138 


0.0454 


7.24 


0.168 


2.55 


73.5 , 221 


M3 


3 


0.147 


0.0417 


6.81 


0.099 


2.15 


76.9 , 231 


M4 


4 


0.150 


0.0361 


6.67 


0.063 


1.80 


76.3 , 229 


M4b 


4 


0.151 


0.0465 


6.63 


0.083 


3.07 


81.3 , 244 


M2b 


2 


0.138 


0.0457 


7.25 


0.166 


2.52 


72.1 , 216 


M2c 


2 


0.138 


0.0456 


7.25 


0.166 


2.52 


74.4 , 223 


TK 


2 


0.135 


0.0180 


7.40 


0.064 


1.68 


149 , 446 



The second column gives the angular periodicity of the spiral or the greatest active sectoral 
harmonic m in the simulation. The next three columns give the pattern speed VLp, growth rate 
7, and the corotation radius Rc of the spiral obtained by fitting it's exponential growth using 
the method in Sellwood & Athanassoula (1986). The next two columns provide the spiral's peak 
amplitude and ^ / (^^)^ ^- And the last column gives the duration of the peak 

amplitude, which we measure between the moment the amplitude reaches the trough and the 
moment preceding the peak at which the amplitude is at the same level as the trough amplitude. 
The second number in the peak duration column is the time in physical units using the adopted 
scaling at the end of section 2.1. 



inner Lindblad resonance lose angular momentum, only a 
tiny fraction end up on retrograde orbits. Note also that 
((AL^)^)^^^ ^ 9.0 X lO'^RiVQ in simulation T, where non- 
axisymmetric forces were eliminated, showing that changes 
in Lz due to orbit integration and noise errors are tiny in 
comparison to those caused by the spiral. 

For each particle, we determine the initial value of 



1 2 1 2 2 



(8) 



where Vz is the vertical velocity component at distance z 
from the midplane and v is the vertical frequency measured 
in the midplane at the particle's initial radius R. For parti- 
cles whose vertical and radial oscillations are small enough 
to satisfy the epicycle approximation, Z — Ez, epi the energy 
of its vertical oscillation; the vertical potential is roughly 
harmonic for \z\ < OARi. Even though the epicycle approxi- 
mation is not satisfied for the majority of particles, we com- 
pute an initial notional vertical amplitude 



(9) 



for them all. The value of defined in this way, i.e. at the 
initial moment only, yields a convenient approximate rank- 
ing of the vertical oscillation amplitudes of the particles, 
although it is clear that in most cases C < Zma^, the maxi- 
mum height a particle may reach. Note that since each disc 
component in our models has an initial thickness that is in- 
dependent of radius, the distribution of C is independent of 
Lz{t = 0). 

Notice from Fig. |5] that changes for the thick disc par- 
ticles are only slightly smaller than those for thin disc par- 
ticles. In order to emphasize this point, the bottom panel 
displays only those thick disc particles for which ^ > 1.2i?i, 
or one scale height of the thick disc, revealing that even for 
these particles changes can be almost as large as those in 
the thin disc. 

Fig. |6] shows the evolution of the radial velocity disper- 
sion for the thin and thick discs. As expected ( Sellwood 
Binney||2002| , the large angular momentum changes near 
corotation cause little heating. Some heating occurs near 



the Lindblad resonances, and is greater near the inner reso- 
nance, again as expected. 



3.2 Distribution of angular momentum change 

Table |3] lists the root mean square, maximum positive, and 
maximum negative changes in angular momentum for both 
all the particles and only those that satisfy > zq in each 
disc. 

Fig. [T] shows the 5th and 95th percentile values of the 
changes in angular momentum as a function of initial no- 
tional vertical amplitude, ^ (^q- [s]), using bin widths of 
O.SRi. The affect of radial migration seems to decrease al- 
most linearly with increasing 

Radial migration should also be weaker for particles 
having larger radial oscillations or epicycles of large am- 
plitude. As reasoned by [Sellwood Binney ^ (2002), par- 
ticles on more eccentric orbits cannot hold station with 
a steadily rotating spiral, because their angular velocities 
vary significantly as they oscillate radially. Fig. [8] which 
is for only those particles in the angular momentum range 
2.5RiVo ^Lz{t = 64) ^ lO.Oi^iVo, confirms that the largest 
angular momentum changes occur among particles having 
the least eccentric orbits. It shows angular momentum 
changes from, and eccentricities at, time t = 64Ri/Vo, which 
is after the model has settled from its mild initial imbalance, 
but before any substantial angular momentum changes have 
occurred. We define eccentricity as e = {Ra — Rp)/ {Ra-\-Rp), 
in which Ra and Rp are respectively the initial apo-centre 
and peri-centre distance of the orbit of a particle having the 
same Lz, but whose motion is confined to the midplane of 
the axisymmetric potential. 



^ The energy cut-off we apply (see after eq. 4) eliminates the 
most eccentric orbits from this plot. 
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Figure 5. Angular momentum changes of particles in M2 as a 
function of their initial angular momenta. The top panel is for 
thin disc particles, the middle for all thick disc particles, and 
the bottom for thick disc particles having notional vertical am- 
plitude C > 1.2i?i. The horizontal red line denotes zero change. 
The vertical lines mark the Lindblad resonances (dashed green) 
and corotation (solid blue). The solid blue line with a slope of 
—2 illustrates the locus of particles whose changes would be sym- 
metric about corotation. The dashed red line of slope —1 shows 
the ALz = —Lz{t = 0) locus below which particles end up on 
retrograde orbits. 
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Figure 6. The radial variations of aR in the thin (bottom curves) 
and thick (top curves) discs of M2 at the same five different times 
as used in the top left and right panels of Fig. [2] with the same 
line colour and style coding. The horizontal dot-dashed curves 
show the theoretical initial dispersions from the untapered discs, 
while the vertical lines mark the principal resonances of the spiral, 
colour and style coded as in Fig. [s] 
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Figure 7. The 5th and 95th percentiles of the angular momentum 
changes in the thick discs of simulations M2 (dot-dashed black), 
M3 (short-dashed red), M4 (solid green), the massive thick disc 
of M4b (long-dashed blue), and the massless thick disc of M4b 
(dotted orange) as a function of notional vertical amplitude C- 
Although we have smoothed the curves, the rising noise with in- 
creasing C is caused by the decrease in the number of particles 
per bin. The curves stop when the number of particles per bin 
drops below twenty. 



© 2011 RAS, MNRAS 000, [Tp^ 



Radial Migration in Thick Discs 9 



> 

N 
< 



1 ' r 



~i ' r 



— ; 




Table 3. Changes in angular momentum resulting from the spi- 
rals and bar. 
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Figure 8. Angular momentum changes of thick disc particles be- 
tween the final time and t = 64i^^/Vb having 2.5RiVo ^ Lz(t = 
64) ^ lO.ORiVo in simulation M2 as a function of their eccentric- 
ities at t = 64:Ri/Vo. The plot for the thin disc is quite similar. 
The red horizontal line shows zero angular momentum change. 



3.3 Effect of radial migration on vertical 
oscillations 

In the previous subsection, we showed how the particles' 
angular momentum changes vary with initial amplitude of 
vertical motion. Here, we present the converse: how the ver- 
tical oscillations are affected by the radial excursions. 

Since the notional amplitude of vertical motion ^ (eq.joj) 
is accurate only in the epicycle approximation, we deter- 
mine a particle's actual maximum vertical excursion, Zma^, 
by integrating its motion in a frozen, azimuthally-averaged 
potential for many radial periods. We do this twice for each 
particle in simulation M2, at time t — QARi/Vo starting from 
the particle's phase space coordinates in the frozen potential 
at that moment and again at the final time t — 387.2i?^/Vb. 

Fig.[9]shows that, on average, the vertical excursions of 
particles increase for those that move radially outwards and 
decrease for those that move inwards. This is expected, be- 
cause restoring forces to the mid-plane are weaker at larger 
radii. For both discs, the mean (solid green) and median 
(dashed blue) curves show a roughly constant Azmax for a 
wide range of AL^. For the thick disc, this constant A^max 
is about twice as great as that for the thin. 

While the large majority of the particles lie in the con- 
toured region, the outliers exhibit significant substructure. 
The dense group of points near a line of slope of —1 in the 
second quadrant are particles with initial home radii within 
and near the inner m — 2 vertical resonanc^ which lies at 
the radius 2.23Ri, not far from the radial inner Lindblad 
resonance at 2.12Ri. Usually u ^ which causes vertical 
resonances to be found significantly farther from corotation 



^ Vertical resonances occur when m(Q — Qp) = zbi^ 
© 2011 RAS, MNRAS 000, [lp2] 



Disc 



R^Vo 



2\l/2 



Max ALz 

R^Vo 



Max -ALz 

RiVo 



M2 


thin 


1.59 


1.56 


4.44 


4.37 


-4.83 


-4.79 




thick 


0.97 


0.85 


4.38 


3.91 


-4.78 


-4.46 


M3 


thin 


0.92 


0.89 


2.79 


2.72 


-2.67 


-2.61 




thick 


0.54 


0.43 


2.76 


2.22 


-2.60 


-2.12 


M4 


thin 


0.62 


0.60 


2.03 


1.97 


-2.11 


-2.02 




thick 


0.37 


0.30 


1.93 


1.65 


-1.99 


-1.59 


M4b 


thin 


0.80 


0.78 


2.25 


2.18 


-2.17 


-2.08 




thick 


0.47 


0.40 


2.21 


1.84 


-2.09 


-1.83 




massless 


0.38 


0.27 


2.17 


1.43 


-2.05 


-1.48 


M2b 


thin 


1.58 


1.54 


4.45 


4.36 


-4.57 


-4.43 




thick 


0.95 


0.84 


4.33 


3.76 


-4.52 


-4.02 




massless 1 


1.08 


1.06 


4.42 


4.32 


-4.43 


-4.43 




massless 2 


1.06 


1.02 


4.41 


4.19 


-4.54 


-4.22 




massless 3 


1.04 


0.98 


4.41 


4.08 


-4.47 


-4.27 




massless 4 


1.02 


0.95 


4.48 


4.11 


-4.43 


-4.15 




massless 5 


0.89 


0.72 


4.35 


3.46 


-4.40 


-3.67 




massless 6 


0.83 


0.64 


4.39 


3.11 


-4.29 


-3.64 




massless 7 


0.79 


0.58 


4.05 


3.06 


-4.30 


-3.37 


M2c 


thin 


1.57 


1.53 


4.46 


4.35 


-4.58 


-4.44 




thick 


0.96 


0.85 


4.33 


3.76 


-4.50 


-4.03 




massless 1 


1.19 


1.06 


4.37 


3.75 


-4.42 


-4.11 




massless 2 


1.08 


0.95 


4.37 


3.81 


-4.46 


-4.01 




massless 3 


U.oO 


U. / D 


A Q'V 
4.0 / 


O.DO 


-4.4y 


A no 

-4.uy 




massless 4 


0.80 


0.72 


4.35 


3.70 


-4.28 


-4.09 




massless 5 


0.75 


0.68 


4.22 


3.79 


-4.33 


-3.89 




massless 6 


0.72 


0.65 


4.18 


3.56 


-4.51 


-3.81 




massless 7 


0.69 


0.63 


4.21 


3.51 


-4.29 


-3.81 


TK 


thick 


0.78 


0.75 


3.37 


3.26 


-2.88 


-2.77 


UC 


thin 


2.65 


2.31 


13.97 


11.60 


-10.43 


, -9.48 




thick 


1.95 


1.75 


12.71 


10.33 


-10.15 


-10.04 


UCBl 


thin 


3.15 


2.76 


16.42 


12.41 


-13.22 


-11.12 




thick 


2.34 


2.10 


16.00 


13.80 


-12.00 


-10.92 


UCB2 


thin 


3.53 


3.46 


19.23 


19.14 


-16.74 


-16.74 




thick 


2.54 


2.28 


19.31 


17.40 


-15.83 


-13.63 



For each disc of every simulation, the first numbers in the third, 
fourth, and fifth columns give the root mean square, maximum 
positive, and maximum negative changes in angular momentum 
respectively. We calculate these for all the particles except those 
that escaped the grid for which the initial angular momentum lies 
in an interval of the spiral's main influence around its corotation. 
This interval, in values in terms of RiVo, is [2.5, 10] for simulations 
M2, M2b, and M2c, [4.5,9.5] for M3, [5.0,8.5] for M4 and M4b, 
and [3.5,9.5] for TK. We do not conflne Lz(t = 0) to such an 
interval for simulations UC, UCBl, and UCB2 since the influence 
of their spirals and bar span almost the entire range. The second 
numbers in the last three columns give the same results but only 
for particles having notional vertical amplitude > zq. 

than the radial resonances, but in our case the inner taper 
reduces the inner surface density so that u ^ k in this part 
of the disc. Thus, these particles are scattered vertically at 
the inner vertical resonance at the same time as they lose 
angular momentum at the radial inner Lindblad resonance. 
Another outlying group can be seen in the first quadrant, 
with large Azmax for small AL^; these particles have quite 
eccentric initial orbits that have particular radial phases just 
before the spiral saturates. Either they are at their pericen- 
tres and lie near corotation just trailing either spiral arm, 
or they are at their apocentres and surround the outer ends 
of the spiral arms. Being at these special locations when 
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Figure 9. Changes in maximum vertical excursions of particles 
in the thin (top) and thick (bottom) discs of M2 versus changes in 
their angular momenta. Contours are linearly spaced in number 
density and we plot only those points that lie outside the lowest 
contour. The mean (solid green) and median (dashed blue) show 
systematic variations, as expected. Note that the vertical scales 
differ in the two plots. 



the spiral wave is strongest, gives them instantaneous angu- 
lar frequencies about the centre that cause them to experi- 
ence more nearly steady, and not oscillatory, torques from 
the perturbation, leading to some angular momentum gain. 
Thus these particles move onto even more eccentric orbits 
and their increased Zmax occurs at their new larger apocen- 
tres. 



3.4 Effects of disc thickness and radial velocity 
dispersion 

Our somewhat surprising finding from simulation M2 is 
that angular momentum changes in the thick disc are only 
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Figure 10. The filled symbols show ^(AL^)^^ as a function 
of vertical thickness (bottom axis) of particle populations of sim- 
ulation M2b. The black squares are for all the particles in the 
angular momentum range 2.5i?iVb ^ Lz{t = ^) ^ ICOi^^Vb, 
while the red triangles are for only those particles having C, > zq. 
The line is a least-squares fit to the black squares of the form 
-2o/5.95i?i rpj^g green plus symbols show the 

\l/2 



variation of ^(AL^)^^ ' with initial radial velocity dispersion 
(top axis) of all the particles in the same Lz{t = 0) range from 
populations in simulation M2c. The blue crosses are for only those 
with Q> ZQ = 1.2Ri. 



slightly smaller than those in the thin, and also in the razor- 
thin disc of ( Sellwood Binney|2002| . Despite the fact that 
thick disc particles both rise to greater z heights and have 
larger epicycles, on average, than do thin disc particles, we 
observe only a mild decrease in their response to spiral forc- 
ing. This finding suggests that the potential variations of a 
spiral having a large spatial scale, such as the m — 2 spiral 
mode in simulation M2, couple well to particles having large 
vertical motions and epicycle sizes. 

To provide more detailed information about how the 
extent of angular momentum changes vary with disc thick- 
ness, we added seven test particle populations to some sim- 
ulations. In M2b, all test particle populations have the same 
initial radial velocity dispersion as the massive thick discs, 
but have scale heights in the range 0.5i?i ^ 2;o ^ 2ARi. Sim- 
ulation, M2c, employs seven test particle populations having 
the same scale height {zo = 1.2Ri) as the massive thick disc, 
but with differing aR. 

Being test particles, they do not affect the dynamics 
of the spiral instability, but merely respond to the poten- 
tial variations that arise from the instability in the massive 
components. These simulations have identically the same 
physical properties as M2 but lower numerical resolution as 
summarized in Table [l] the reduced numerical resolution re- 
mains adequate since the fitted spiral mode is little changed 
from that in simulation M2 (Table [2|. 
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Variations of (^{ALz)^Y^'^ with both disc thickness, at 
fixed radial velocity dispersion, and of radial dispersion at 
fixed thickness, are displayed in Fig. [lO] The decrease is 
somewhat more rapid in subpopulations of particles that 
start with a notional vertical amplitude C > ^o, as seems 
reasonable. The fitted line indicates that (^(ALz)'^^^^'^ de- 
cays approximately exponentially with disc thickness with a 
scale that can be related to theory - see §4.3| Table |3] gives 
both the plotted root mean square, as well as the maximum 
positive, and maximum negative ALz for all the particles 
of each population and separately for those with ^ > ^2:0 of 
each disc. 

We caution that the information in Table |3] and in 
Fig. [To] quantifies how the responsiveness of a test parti- 
cle population scales when subject to a fixed perturbation. 
Self-consistent spiral perturbations may differ in strength, 
spatial scale, and/or time dependence causing quite differ- 
ent angular momentum changes. 

3.5 Thick disc only 

For completeness, we also present simulation TK, which has 
a single, half-mass active disc with a substantial thickness. 
We inserted the same initial groove, and restricted forces to 
m = 2 only. 

It is interesting that a groove in the thick disc still cre- 
ates a spiral instability, but one that grows less rapidly and 
saturates at a lower amplitude than we found in M2. As a 

consequence, the spread in (^{ALz)'^^^'^ about 1/2 as large 
as in M2. Again we find that the radial migration is reduced 
by increasing disc thickness, but not inhibited entirely. 



to local theory, depends on the vigour of the swing amplifier 
( Toomre|1981| . Thus to generate a similar spiral disturbance 
with m > 2, we need to hold the key parameters X and Q 
at similar values. For an m-armed disturbance in a thin, sin- 
gle component Mestel disc with active mass fraction /, the 
locally-defined parameter 

(11) 

fm 

is independent of radius. We therefore scale the active mass 
fractions in both components as / oc - recall that we 
used a half-mass disc for m = 2. In order to preserve the 
same Q value (eq. |3|, the radial velocity dispersion of the 
particles also has to be reduced as an (x f. 

Simulations M3 and M4 therefore have lower surface 
densities and smaller velocity dispersions in order to support 
similarly growing spiral modes of sectoral harmonic m — 3 
and 4 respectively. A further simulation M4b is described 
below. These models are all seeded with a groove of the 
same form (eq. [g]) and parameters as that in M2. 

Note that the instability typically extends between the 
Lindblad resonances, which move closer to corotation as m 
is increased. In the Mestel disc, the radial extent of the mode 
varies as 

i^OLR _ m + V2 ^^2^ 
RiLR m — a/2 

i.e. i^oLR/i?iLR ^ 5.8, 2.8 k 2.1, for m = 2, 3 & 4 respec- 
tively. Since we have also decreased the in-plane random 
motion in proportion to the surface density decrease, the 
decreased size of the in-plane epicycles somewhat compen- 
sates for the smaller scale of the mode, although the ratio is 
not exactly preserved. 



4 OTHER SIMULATIONS 

The potential of a plane wave disturbance in a thin sheet 
having a sinusoidal variation of surface density Sa in the 
x-direction is 

$,(x,z) = -?^e^^-e-l^^l (10) 
1^1 

where k is the wavenumber (BT08, eq. 5.161). The exponen- 
tial decay away from the disc plane is steeper for waves of 
smaller spatial scale, i.e. larger While spirals are not sim- 
ple plane waves in a razor-thin sheet, this formula suggests 
that we should expect radial migration to be weakened more 
by disc thickness for spirals of smaller spatial scale or higher 
angular periodicity - note that the azimuthal wavenumber, 
k(f) — m/R. The simulations in this section study the effect 
of changing this parameter. 

4.1 Spirals of different angular periodicities 

We wish to create spiral disturbances that are similar to 
that in M2, but have higher angular periodicities. Since we 
expect migration to depend not only on the spatial scale 
relative to the disc thickness, but also the peak amplitude, 
Sa, and perhaps also the time dependence, we try to keep 
as many factors unchanged as possible. 

The mechanics of the instability seeded by a groove de- 
pends strongly on the supporting response of the surround- 
ing disc ( jSellwood h Kahn||199l| , which in turn, according 



4.2 Results for m > 2 

Table [2] gives our estimates of the spiral properties in each 
simulation; uncertainties in the measured frequencies are 
typically < 2% (Se llwood k AthanassouIa ''1986). The pat- 
tern speeds of these instabilities do not change much with 
m, except that we find the radius of corotation lies closer 
to the groove centre L*/Vo = 6.5i?2, as expected. The table 
also gives the time during which the amplitude of the wave 
is equal to or greater than that at the post-peak trough, 
which again does not vary much with angular periodicity. 
However, both the growth rates and the peak amplitudes of 
the modes decrease from M2, to M3 and M4. 

Table [3] includes the root mean square, maximum posi- 
tive, and maximum negative angular momentum changes for 
M3, M4, and M4b, measured in each case to the moment at 
which Am passes through the first minimum after the mode 
has saturated. Comparison with M2 reveals that increasing 
m causes a roughly proportionate decrease in the angular 
momentum changes, in part because the saturation ampli- 
tude is lower, but perhaps also because the spatial scale is 
reduced. Note that again some thick disc particles in both 
M3 and M4 have ALz values almost as large as the greatest 
in the thin disc, as was also the case for M2. 

Fig.|7|shows 5th and 95th percentile values of ALz ver- 
sus initial notional vertical amplitude C, for simulations M3, 
M4, and M4b. Compared to the curve of M2, the ALz values 
are smaller for greater m - i.e. the extent of radial mixing is 
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substantially lessened. Note that this difference could have 
a variety of causes, such as the lower limiting amplitude of 
the mode, or possibly the different growth rate of the mode, 
and/or the different disc thickness relative to the spatial 
scale of the mode. 

This last factor is one we are able to eliminate. The 
thickness of the discs of simulations M2, M3, and M4 were 
held fixed as we increased m and reduced the surface density. 
In order to eliminate a change in the ratio of disc thickness 
to spatial scale of the mode, we ran a further simulation 
M4b with the same in-plane parameters as in M4, but with 
half the disc thickness. We also halved the gravity softening 
length and the vertical spacing of the grid planes. 

This change restores the growth-rate of the mode in run 
M4b to a value quite comparable to that in simulation M2 
(Table [2]). The saturation amplitude, while larger than in 
simulation M4, is still about half that in M2. In addition, 
we added a test particle population to simulation M4b that 
has parameters identical to those of the massive thick disc, 
including an or q, except its vertical scale height zq is that 
of the thick disc of M2, M3, and M4. Again as expected. 
Table [3] reveals that (^{ALz)'^^^'^ is substantially lower in 
the massless disc than in the thinner, massive disc. 



4.3 Comparison with theory 

In order to make sense of these results, we here compare 
with the theoretical picture developed by |Sellwood Binney| 
(20021. 



First, we eliminate the possibility that the spiral in the 
simulations with higher m is "on" for too long for optimal 



migration. Sellwood & Binney (2002) argue that efficient 



mixing by the spiral requires the duration of the peak am- 
plitude be less than half the period of a horseshoe orbit, so 
that each particle experiences only a single scattering. They 
show that the minimum period of a horseshoe orbit varies as 
1^0 1 where the potential amplitude of the spiral pertur- 
bation at corotation varies with the spiral density amplitude, 
Sa, and sectoral harmonic as |^o| oc T^a/m (eq. 10). Thus 



the weaker density amplitude that we find with higher m 
implies that the minimum periods of the horseshoe orbits 
are greater, and the condition for efficient mixing is more 
strongly fulfilled for m > 2. 

Theory also suggests that an m-dependence of the peak 
amplitude is unavoidable if the spiral saturates due to the 
onset of many horseshoe orbits, as was argued in |Sell-| 



wood & Binney (2002). In their notation, orbits librate 

< p^, where the frequency 
10), we see that 



i.e. are horseshoes - when En 



oc ml^ol"*"^^. Since |^o| oc T^a/ra (eq. 



p oc mHa, suggesting that as m increases, horeshoe orbits 
become important at a lower peak density. Comparing M2 
with M4b, we find (Table [2| the relative limiting amplitudes 
Ea oc m~^, which is consistent with the idea that the spiral 
instability saturates when the importance of horseshoe or- 
bits reaches very nearly the same level. The fixed thickness 
and softening length used in M2, M3 and M4 dispropor- 
tionately weakens the potential of the spiral as m rises, and 
spoils this exact scaling. 

Note also that both (^(ALz)'^Y^'^ and the extreme val- 
ues measured from M4b are almost exactly half those in M2, 
which is also consistent with the horseshoe orbit theory de- 



veloped by Sellwood k Binney (2002). Since they showed 
that the maximum ALz oc |^o| ^ , the relations in the pre- 
vious paragraph require AL^^max oc as we observe. As 
the rms value also scales in the same way, it would seem the 
entire distribution of ALz scales with m in the same way. 
Again, the constant disc thickness prevents this prediction 
from working perfectly for M3 and M4. 

We stress that this scaling holds because we took some 
care to ensure the key dynamical properties of the disc were 
adjusted appropriately. Spirals in a disc having a different 
responsiveness would have both a different growth rate and 
probably also peak amplitude, and the behaviour would not 
have manifested a simple m-dependence. 

Finally, we motivate the exponential fit to the variation 
of {{ALzf) with zo, for the same m — 2 spiral distur- 
bance. We have already shown that ^(AL^)^^"^^^ oc |^o|"^^^, 
and have argued that the spiral potential decays away from 
the mid-plane as e"''^^' (eq. 10), with /c^ = m/R. Naively, we 
could set z = 2^0, k = kci)=m/R, m = 2, and R = Rc, since 
angular momentum changes are centred on corotation, lead- 



ing to {{ALzf) oc e-^o/^^ The fitted scale is 5.95i?^, 
which is somewhat smaller than Rc = 7.24i?i. The domi- 
nant cause of this discrepancy is probably that spiral is not 
a plane wave, curvature is important for m = 2, and that its 
wavenumber is larger than /c^ because the spiral ridges are 
inclined to the radial direction; we should therefore expect 
the potential to decay away from the mid-plane rather more 
rapidly, in the sense that we measure. 

We conclude that angular momentum changes scale 
with spiral amplitude in the manner predicted in |Sellwood| 
|fc Binne y] (|2002 ) and, furthermore, the limiting amplitude 
itself is determined by their theory. The variation with disc 
thickness is also in the sense expected from the theory, but 
the quantitative prediction is not exact. 



5 UNCONSTRAINED SIMULATIONS 

Having studied at length the effects of a single spiral wave, 
we now wish to illustrate the effects of multiple spirals. We 
present three simulations, UC, UCBl, and UCB2, that have 
no initial groove and the initial positions of the particles 
are random (i.e. not a quiet start) since we wish spirals 
to develop quickly from random fluctuations. We include 
gravitational disturbance forces from all sectoral harmonics 
^ m ^ 8, except m = 1, which we omit in order to avoid 
imbalanced forces from a possibly asymmetric distribution 
of particles in a rigid halo with a fixed centre. 

All the many simulations of this type that we have run 
have ultimately developed a strong bar at the centre. We 
here compare radial migration in three separate cases: sim- 
ulation UC avoids a bar for a long period and supports many 
transient spirals, while simulations UCBl and UCB2 have 
identical numerical parameters but form a bar quite early. 
In all three cases, the combined thin and thick discs have 
a smaller active mass fraction, / = 0.44, compared with 
/ = 0.55 for simulation M2. A smaller active mass and a 
larger gravitational softening length help to delay bar for- 
mation but also weaken the m — 2 spiral amplitudes some- 
what. The numerical parameters of these simulations are 
also listed in Table [T] Since bar formation in these models 
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Figure 11. The time evolution of A2/A0 in simulations UC 
(black), UCBl (red), and UCB2 (green). The top axis shows time 
scaled to physical units using the adopted scaling. 



is stochastic {e.g. Sellwood & Debattista 2009'), the differ- 
ent bar-formation times simply arise from different initial 
random seeds. 



5.1 Multiple spirals only 



The black curve in Fig. 11 shows the evolution of A2/A0 
in simulation UC. The spiral amplitudes rise to the point 
at which significant particle scattering begins by t ^ 
1500i?i/Vb and we present the behaviour up to time 
3 500i?i/Vb shortly before a bar begins to form. The pe- 
riod IbOORi/Vo < t < 3 500Ri/Vo during which particles 
are scattered by spiral activity corresponds to ~ 6.0 Gyr 
with our adopted scaling. 

Fig. 1 12| illustrates the m — 2 power spectrum of distur- 
bances. Each horizontally extended peak indicates a spiral 
of a particular pattern speed. Some 20 transient spirals of 
significant amplitude occur during this period spanning a 
wide range of angular frequencies and corotation radii. 

Since there are numerous disturbances with well scat- 
tered Lindblad resonances, random motion rises generally 
over the disc, as reported in Fig. [13] in contrast to Fig. [6] 
which shows the localized heating at the ILR of the single 
spiral case. 

The larger changes in Lz than those for the single spiral 
case are evident from the dot-dashed black curves of Fig.[l4] 
which show the 5th and 95th percentiles of AL^^ versus initial 
notional vertical amplitude C- Although the shapes of these 
curves are similar to those in Fig. |7| they are asymmetric 
about zero indicating that gains are larger than losses. 

The angular momentum changes illustrated in Fig. [15] 
reflect mostly the effects of the latest spirals that devel- 
oped in the simulation in each Lz{t — 0) region. Contours 
of changes to the distribution of home radii for the particles 



(dotted black and solid green in Fig. 16 ) reveal that particles 
from all initial radii can move to new home radii, but that 
the changes are greatest around the mid range of initial radii 
where corotation resonances are more likely. Some 0.003% 
of the particles in the thin disc and 0.03% in the thick end 
up on retrograde orbits. 




R/R, 



Figure 12. The power spectrum of m = 2 density variations in 
simulation UC. The solid curve indicates the radius of corotation 
for the given frequency, while the dashed curves show the radii of 
the Lindblad resonances. 
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Figure 13. Radial variations of ctr in simulation UC for the thin 
(solid) and thick (dashed) discs. The curves are drawn at equal 
time intervals and (Jr rises monotonically in the inner and outer 
discs. 



Taken together with the results from the single spiral 
simulations described above, we again find that the changes 
in angular momenta and home radii are almost as great in 
the thick disc as in the thin. That is, radial migration is 
weakened only slightly by disc thickness. 

The top row of Fig. [17] shows changes in the UC par- 
ticles' maximum vertical excursions versus changes in their 
angular momenta. We find that, on average, Zmax increases 
except for the greatest losses in angular momentum. This is 
different from the single spiral case, in which the contours 
and the mean and median curves are centred on the origin 
(Fig. [9]). This overall extra increase in vertical amplitude 
probably comes from the net heating by vertical resonances 
that the multiple transient spirals induce by the end of the 
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Figure 14. Same as Fig.|7|for simulations UC (dot-dashed black) 
and UCBl (long-dashed red). The solid black and short-dashed 
red curves show the average angular momentum changes for UC 
and UCBl respectively, and the dotted black and medium-dashed 
red ones show the median changes. 



simulation. Nevertheless, the trend of the mean and median 
with ALz remains as in M2. An exception is again the group 
of points along the slope of —1 in the second quadrant. It 
is much more pronounced for the thin disc than in M2. The 
particles contributing to this feature are still initially from 
the inner region of the disc, but this region is more extended 
in UC since the inner vertical resonances occur at various 
radii for the numerous transient spirals. As in M2, changes in 
^max remain about twice as great for the thick disc particles 
as for those of the thin for the same ALz. 

Although the scale height of outward migrating parti- 
cles does increase somewhat, as expected, the changes are 
not substantial enough to cause the thickness of outward 
migrating thin-disc particles to approach the scale height of 
the thick disk. The changes in the vertical motion of inward 
migrating particles are more minor than those for outward 
migrators. Thus we do not observe much of a tendency in our 
models for evolution to cause a significant degree of blurring 
between the separate populations. 



5.2 Multiple spirals with a bar 

We here study radial migration in two simulations, UCBl 
and UCB 2, that formed bars at an early stage of their evo- 
lution and compare them with simulation UC that did not 
form a bar for a long period. As noted in the introduction, 
it has long been known that bar formation causes some of 
the largest changes to the distribution of angular momen- 
tum within a disc. Here our focus is on the consequences of 
continued transient spiral behaviour in the outer disc long 
after the bar formed, which has not previously received much 
attention, as far as we are aware. 
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Figure 15. Same as Fig.[5]for simulation UC. 
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Figure 16. The distributions of final home radii versus initial 
home radii for all the particles in simulation UC. The particle 
density in this plane is estimated using an adaptive kernel, the 
dotted black and solid green contours represent the thin and thick 
discs respectively, and the red line shows zero change in i?^ . Con- 
tour levels are chosen every 10% of the thick disc's maximum 
value from 5% to 95% for both thin and thick discs. The arrow 
light blue and dashed orange contours represent the same for the 
thin and thick discs of UCBl respectively. The region inside the 
bar {Rhit = 3500) < lOi?^) is omitted. 



The red curve of Fig. 11 shows the evolution of A2/AQ 
in UCBl. The spirals become significant around the time 
~ 1000i?^/Vb, the bar forms around ~ 2 000i?i/Vb, and we 
stop the simulation at the same final time t — 3 500i?^/Vb 
as UC. Thus, the time interval of significant scattering is 
roughly 2 500i?i/Vb in length with a bar being present for 
the last ^ 1500i?i/Vb, which correspond to 7.5 Gyr and 
4.5 Gyr respectively. With our suggested scaling, the bar 
has a pattern speed of 31.8 km s~^ kpc~^ and corotation 
Rc ^ 7.5 kpc. (This scaling makes the bar substantially 
larger than that in the Milky Way.) 

The long-dashed red curves in Fig. ^] show that the 
bar enhances the changes in Lz somewhat over those that 
arise due to spirals alone. Note that we measure the in- 
stantaneous value of Lz of each particle and it should be 
borne in mind that it changes continuously in the strongly 
non-axisymmetric potential of this model, especially so for 
particles in or near the bar. 

For this reason, we cannot extend the light blue (marked 



with arrow heads) and dashed orange contours in Fig. 16 
to small home radii at the later time in the strongly non- 
axisymmetric potential of the bar. However, a clear asym- 
metry can be seen; the distribution is biased above the zero 
change red line, indicating a systematic outward migration 
in the outer disc. A similar asymmetry can be seen without 
a bar in UC (dotted black and solid green contours), but the 
formation of the bar makes it more pronounced. Since total 
angular momentum is conserved in these simulations, there 
is a corresponding inward migration in the inner regions (see 
Fig.|5f. 

Aside from the bar region, where systematic non- 
circular streaming biases the rms radial velocities. Fig. [18] 
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Figure 17. Same as Fig.jojfor simulations UC (top row), UCBl 
(middle row), and UCB2 (bottom row). The left column shows 
A^^max as a function of ALz for the thin discs and the right for 
the thick. 
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Figure 18. Same as Fig.flSlfor UCBl. 



shows that the bar does not appear to cause much extra 
heating over that in the unbarred case. 

Fig. |19| again illustrates that changes in Lz in the outer 
disc are more substantial in this barred model than in the 
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unbarred case (Fig. 15), but the larger spread in ALz in the 
inner disc is partly an artifact of using the instantaneous 
values of Lz at later times. Because we use the instantaneous 
Lz in the barred potential in this Figure, a particle in the 
region where ALz < —Lz need not necessarily have been 
changed to a fully retrograde orbit. 

The asymmetries in Figs. ^] & ^] about the red lines 
of zero change are not caused by bar-formation alone, as 
the distributions (not shown) are approximately symmetric 
immediately after this event. Rather, we find they appear to 
be caused by multiple migrations of the same particles by 
separate events. 

Our second barred simulation, UCB2, again differs from 
UCBl and UC only by the initial random seed. It formed 
a stronger bar, but a little later than in UCBl, as shown 
by the green curve of Fig. |11| Significant scattering occurs 
for the last ^ 1800i?i/Vo time units (5.4 Gyr) of which 
the bar is present for the last ^ 1 lOOi^i/Vb (3.3 Gyr). The 
pattern speed (33.6 km s~^ kpc~^) and corotation radius 
7.3 kpc) for our adopted scaling, are similar to the values 
in UCBl. 

We find the extent of radial migration (bottom two rows 
of Table [3| is further increased by the stronger bar, but 



10.2 



19] & [16] (not shown) are 
but the 



not by much. Variants of Figs, 
qualitatively similar with slightly larger changes, 
outer disc is again dominated by scattering due to the latest 
few spirals. 

In both these models, therefore, we see that the forma- 
tion of a bar does indeed increase the net angular momentum 
changes. However the overall behaviour is similar to that of 
scattering by transient spirals without a bar. This is differ- 
ent from the effect found by Brunetti et al. (2011)), in which 
essentially all the angular momentum changes occurred dur- 
ing bar formation. 

The lower two rows of Fig. [iTj show Azmax versus ALz 
for UCBl and UCB2 respectively. The presence of a bar 
yields a much denser and more extended feature in the sec- 
ond quadrant, which is also visible in for thick disc. For 
UCB2, its contribution is so great that the mean ZA^^max 
keeps increasing with greater angular momentum loss. This 
extra apparent vertical heating is probably caused by the 
buckling of the bar. Unlike for M2, UC, and UCBl, we find 
that for positive ALz in UCB2, the mean and median curves 
keep rising for larger ALz- 



6 A CONSERVED QUANTITY? 

Radial migration results from angular momentum changes 
near corotation, which we have shown to be somewhat 
weakened by increased vertical motion. [Schonrich Bin^ 



ney 



(2009a), in their model of radial mixing in thin and 
thick discs, assumed that vertical and radial motions are 
decoupled and that vertical energy is conserved as stars mi- 
grate radially. We here try to identify a conserved quantity 
that can be used to predict vertical motion when particles 
suffer large changes in Lz, by comparing various measures 
of vertical amplitude at the initial and final times in sim- 
ulations with a single spiral. Note that the "initial" value 
of each quantity we discuss in this section is measured at 
t — 64i?i/Vb, which avoids possible effects of the settling of 
the model from its mild initial imbalance. Also, all quantities 
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Figure 19. Same as Fig. |15| for simulation UCBl. The vertical 
blue line marks the approximate corotation resonance of the bar. 
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Figure 20. The top left panel shows a typical thick disc orbit in 
the meridional plane, and the other three panels show how Ez^r 
varies with time (top right), radius (bottom left), and z-height 
(bottom right). The particle, which has a home radius of 6.5Ri, 
a total energy of —2.0Vq, a midplane eccentricity of 0.15, and 
a vertical excursion of 2.1Ri, is integrated in the frozen initial 
axisymmetric potential of simulation M2. 



are computed in an azimuthally averaged potential, to elim- 
inate variations with spiral phase, which remain significant 
at the final time. 

We focus on particles whose initial home radii lie in 
an annulus of width G.ORi centred at corotation of the spi- 
ral in M2, and measure quantities for only one particle per 
quiet start ring, meaning one in every twelve. This results 
in ^ 40 000 particles from the thin disc and ^ 73 000 from 
the thick. Although the entire thick disc is represented by 
just 50% more particles than is the thin, the effects of 
the inner and outer tapers, together with the groove in 
the thin disc cause the number of particles in the range 
4.0Ri ^Rh{t = 64) ^ lO.ORi to be some 82% larger for the 
thick disc than for the thin. While we endeavour to measure 
each quantity in this section for the same set of particles, 
from both simulations M2 and T, we have been able to es- 
timate some quantities for only a subset of these particles, 
as noted below. 



6.1 Vertical energies 

We have tried a number of different way to estimate the 
energy of vertical motion. A simple definition might be 



E.,ii= -vl + <i>iR,z)-<S>iR,0), 



(13) 



with R being the instantaneous radius of the particle. How- 
ever, Fig. [20| shows that Ez,r defined this way varies by some 
30% as a particle oscillates in both radius and vertically in 
a static axisymmetric potential. The orbit of this particle 
is multiply periodic and clearly respects three integrals (as 
we will confirm below), in common with many orbits in ax- 
isymmetric potentials (BT08, section 3.2). Although such 
orbits can be described by action-angle variables, which im- 



ply three decoupled oscillations, the energy of vertical mo- 
tion is clearly not decoupled from that of the radial part of 
the motion. In particular, the bottom left panel shows that 
Ez,R varies systematically with radius, which refiects in part 
the weakening of the vertical restoring force with the out- 
wardly declining surface density of the disc. This behaviour 
considerably complicates our attempts to compare the ver- 
tical energies at two different times in the same simulation. 

The epicycle approximation (BT08, p. 164) holds for 
stars whose orbits depart only slightly from circular motion 
in the midplane. In this approximation, when a star pursues 
a near-circular orbit near the mid-plane of an axisymmet- 
ric potential, the vertical and radial parts of the motion are 
separate, decoupled oscillations, and the vertical energy is 
Ez,epi (= Z eq.[8| is constant. However, the epicycle approx- 
imation is a poor description of the motion of most particles, 
for which the radial and vertical oscillations are neither har- 
monic nor are the energies of the two oscillations decoupled. 

Table [I] gives the rms values of [^(tfinai) — 
^(tinitiai)]/^(tinitiai) for particlcs in both the thin and the 
thick discs of simulations M2 and T^ Here Y represents 
one of several possible vertical integrals. The first rows 
give the fractional changes in the epicycyclic approxima- 
tion, Y = Ez,epi, that are substantial. Furthermore, they 
are almost as large in simulation T, which was constrained 
to remain axisymmetric, as those in M2 in which substantial 
radial migration occurred, suggesting that the changes are 
mostly caused by the inadequacy of the approximation. 

The second rows in Table H show the rms fractional 
changes in the instantaneous values of Ez,r. While these 
values are slightly smaller than for the epicyclic estimate, 
they are again large for the thin disc and still greater for 
the thick. This is hardly surprising, as the particles in the 
simulation have random orbit phases at the two measured 
times. 

Particles on eccentric orbits generally spend more time 
at radii R > Rh than inside this radius, which biases the in- 
stantaneous measure to a lower value, as shown in the third 
panel of Fig. |20| To eliminate this bias, and to reduce the 



random variations, we evaluate Ez,r^ from eq. (13) at the 
home radius of each particle, which requires us to integrate 
the motion of each particle in the frozen potential of the 
appropriate time until the particle reaches its home radius. 
A tiny fraction, ~ 500 of the ^ 113 000 sample particles, 
never cross R = Rh at the initial time and ^ 200 more 
at the final time; these particles rise to large heights above 
the mid-plane and have meridional orbits resembling that 
shown in the top right of Fig. |A2| but appear to be con- 
fined to R > Rh. The fact that this is possible seems consis- 
tent with the description of vertical oscillations developed 
by [Schonrich Binney| ( |2012| ). The third rows of Table [I] 
give the fractional rms changes in Ez,r^ for all the remaining 
particles; the changes in simulation M2 are a little smaller 
than those of the instantaneous values and significantly so 
in simulation T. 



^ For all quantities in this table, we use the biweight estimator of 
the standard deviation ( [Beers et a/.||1990| , which ignores heavy 
tails. We also discard a few particles with initial values < 10~^ 
in order to avoid excessive amplification of errors caused by small 
denominators. 



© 2011 RAS, MNRAS 000,[lp2l 



18 M. Solway, J. A. Sellwood & R. Schdnrich 



Table 4. Biweight standard deviation of fractional changes in 
various estimates of vertical energy and action. 
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The biweight estimated standard deviation of the fractional 
change in variable Y listed in the third column for both the thin 
and thick discs in simulations M2 and T. We give two values 
for the standard deviation in some rows: the value in the fourth 
column is measured from all the particles, that in the fifth col- 
umn is from only the 48% (in M2) of particles for which AJz is 
calculable. 



Since Ez,r is multiply periodic (Fig. 20), we have also 
estimated an orbit- averaged value, {Ez,r), by integrating 
the motion for many periods until the time average changed 
by < 0.1% when the integration is extended for an addi- 
tional radial period. The fourth rows of Table |4] give the 
rms changes in {Ez,r), which are still considerable. Note 
that we did not compute this time-consuming estimate for 
all the particles, but for only the subset used for other values 
in the fifth column of this table. The choice of this subset is 
described below. 

Generally, we find that none of these estimates of the 
vertical energy is even approximately conserved, except for 
the orbit- averaged energy when the simulation is constrained 
to remain axisymmetric (simulation T). In this case, the po- 
tential at the two times differs slightly due to radial varia- 
tions in the mass distribution - values of this orbit- averaged 
quantity in an unchanged potential would be independent 
of the moment at which the integration begins. 

In all cases, changes are larger in simulation M2 where 
significant radial migration occurs. Fig. [2l] illustrates that 
the changes in the estimated vertical energy correlate with 
ALz. The excesses of particles in the second and fourth 
quadrants indicate that changes have a tendency to be neg- 
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Figure 21. Changes in Ez^epi (first row), Ez,r (second row), 
Ez,Rh (third row), and {Ez^r) (fourth row) versus ALz for the 
thin (left column) and thick (right column) discs of simulation 
M2. Just as in Fig. [o] the horizontal and vertical red lines show 
zero changes, the linearly spaced magenta contours show number 
density, and the bold solid green and dashed blue curves show 
the mean and median changes in each ordinate respectively. Note 
that the vertical scales differ in each plot. 



ative for outwards migrating particles and positive for in- 
wards migrating particles. The significance of this trend is 
discussed below. 



6.2 Vertical actions 

The various actions of a regular orbit in a steady poten- 
tial are defined to be (27r)~^ times the appropriate cross- 
sectional area of the orbit torus (BT08, pp. 211-215). One 
advantage of actions is that they are the conserved quanti- 
ties of an orbit when the potential changes slowly - i.e. they 
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Figure 22. Same as Fig. [21] for changes in Jz,epi (first row), 
Jz\r,cI) (second row), and Jz (third row) of the particles in M2. 



are adiabatic invariants under conditions that are defined 
more carefully in BT08 (pp. 237-238). 
The vertical action is defined as 



Jz = ^ (p zdz, 



(14) 



where z and z are measured in some suitable plane that 
intersects the orbit torus. Since Lz (= Jcf)) is conserved in 
a steady axisymmetric potential, the orbit can be followed 
in the meridional plane (BT08, pl59) in which it oscillates 
both radially and vertically, as illustrated in the first panel 
of Fig. |2Q[ To estimate J^, we need to integrate the orbit 
of a particle in a simulation from its current position in the 
frozen, azimuthally averaged potential at that instant and 
construct the (z, z) surface of section (SoS) as the particle 
crosses Ru with R > 0, say. The integral in eq. ([14| is the 
area enclosed by an invariant curve in this plane. 

Before embarking on this elaborate procedure, we con- 
sider two possible approximations. When the epicycle ap- 
proximation holds, the vertical action is Jz,epi = Ez,epi/i^ 
(BT08 p. 232). The fifth rows of Table [2] for each disc show 
that changes in this quantity are large and again they are 
similar to those in simulation T, confirming once again that 
the epicycle approximation is inadequate. 

Since most orbits reach beyond the harmonic region 
(\z\ < OARi) of the vertical potential, an improved local 
approximation is to calculate 



at the particle's fixed position in the disc. This local esti- 
mate still ignores the particle's radial motion, but gives a 
useful estimate for the average vertical action that is also 
used in analytic disc modeling ( |Binney|2010| . The function 
z{z) at this fixed point is simply determined by the vertical 
variation of $(i?, 0, z), and the area is easily found. We eval- 
uate Jz\r,cI) at the particle's instantaneous position at the 
initial and final times using the corresponding azimuthally 
averaged frozen potential. The rms variation of AJz\R,(f) is 
given in the sixth rows of Table ^ for each disc. We find that 
AJz\R,ct) in the thin discs of both M2 and T are small, sug- 
gesting that this estimate of vertical action is more nearly 
conserved. However, the changes for the thick discs are still 
large, and we conclude that this local estimate is still too 
approximate. 



We therefore turn to an exact evaluation of eq. ( 14) us- 
ing the procedure described in the Appendix. Unfortunately, 
we can evaluate AJz only if the consequents in the SoS form 
a simple invariant curve that allows Jz to be estimated at 
both times. We find that only 48% of the - 113 000 parti- 
cles have closed, concave invariant curves at the initial time 
and only 73% of those retain these properties at the final 
time. Although the thick disc contains more particles than 
the thin, we are able to calculate AJz for a smaller fraction: 
we succeed with ^ 26 000 in the thin disc but only ^ 13 000 
in the thick. 

Column five of Table ^ gives the rms changes of all esti- 
mates of vertical energy and action for only these particles. 
While the fractional rms values of the different energy es- 
timates remain large, they are generally smaller than those 
for the all particles given in the fourth column, but only 
slightly so for Jz\R,(f)' Thus the orbits for which AJz can be 
computed are not a random subset, but are biased to those 
for which energy changes are smaller on average. 

The seventh (final) rows of Table [4] for each disc list the 
rms values of the fractional change in AJz. They are just a 
few percent for particles in simulation T, in which the disc 
is constrained to remain axisymmetric, and are smaller than 
for any other tabulated quantity in the discs perturbed by 
a spiral. The small scatter about the red line of unit slope 
in Fig. [23] indicates that differences in the final and initial 
values of Jz in both discs for M2 are indeed small. 

Fig. [22] shows that changes in all three estimators of 
the vertical actions are uncorrelated with ALz, and that 
the mean and median changes are close to zero. Comparison 
with the systematic changes in Fig. |21| suggest that Jz is, in 
fact, conserved. To see this, note that changes in Lz shift the 
home radius of each particle to a region where the average 
vertical restoring force differs, the amplitudes of the verti- 
cal oscillations increase due to the weaker average restoring 
force when ALz > 0, and conversely decrease for ALz > 0. 
Thus if an estimator of vertical motion is not conserved, 
we should expect a systematic variation with ALz, as we 
observe for all the energy estimators in Fig. |21| The fact 
that there is no systematic variation of the median or mean 
change in vertical action suggest that it is conserved. 

If Jz is a conserved quantity, it may seem puzzling that 
the rms changes in our measurements are not smaller. Both 
Jz,epi and Jz\r,cI) are approximate, which could be respon- 
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Jz(ti)/RiVo 

Figure 23. Comparison of Jz for ~ 40 000 particles from simula- 
tion M2 at the initial and final times. Values are calculated from 
eq. ( |14| . The top panel is on a linear scale, the bottom on a log 
scale to reveal the behaviour for very small actions. The red line 
has unit slope. Thin disc particles are marked in black, thick disc 
in green. 



sible for apparent substantial changes, but changes in the 
"exact" estimate of Jz remain significant, even for simula- 
tion T. It is possible that our estimate of AJz is inaccurate. 
For example, we must eliminate non-axisymmetric structure 
from the potential to compute J^, introducing small errors 
when the potential is mildly non-axisymmetric. Larger dif- 
ferences could arise, however, if the invariant curve changes 
its character, especially since entering or leaving a trapped 
area in phase space is not an adiabatic change. It is possible 
an orbit that appears to be unaffected by resonant islands in 
phase space at the initial and final times could have experi- 
enced trapping about resonant islands for some of the inter- 



mediate evolution, allowing Jz to change even when changes 
to the potential are indeed slow. 

We have verified that Jz is conserved to a part in lO'^ 
in a further simulation with a frozen, axisymmetric poten- 
tial. Note this is not a totally trivial test, since we use the 
A^-body integrator and the grid-determined accelerations to 
advance the motion for ^ 350 dynamical times before mak- 
ing our second estimate. Therefore the non-zero changes in 
simulation T are indeed caused by changes to the poten- 
tial. Even though the potential variations are small, sub- 
stantial action changes indicate that changes to the orbit 
were non-adiabatic, which can happen for the reason given 
in the previous paragraph. Our "initial" Jz values are es- 
timated at t = 64i?i/Vb, when we believe the model has 
relaxed from the initial set up; we therefore suspect that 
the small potential variations are more probably driven by 
particle noise, which can be enhanced by collective modes 



( Weinber'g]|l998| . This hypothesis is supported by our find- 
ing variations in Jz that were about ten times greater than 
in simulation T when we employed 100 times fewer particles. 



7 SUMMARY 

We have presented a quantitative study of the extent of ra- 
dial migration in both thin and thick discs in response to a 
single spiral wave. We find angular momentum changes in 
the thick disc are generally smaller than those in the thin, 
although the tail to the largest changes in each population 
is almost equally extensive. 

We have introduced populations of test particles into 
a number of our simulations in order to determine how 
changes in Lz vary with disc thickness and with radial veloc- 
ity dispersion when subject to the same spiral wave, finding 

.1/^ ....r.-.^Q 



, as shown in Fig. 



an exponential decrease in ^(AL^)^^ 

We find that spirals of smaller spatial scale cause 
smaller changes. When we were careful to change all the 
properties of the model by the appropriate factors, we were 

able to account for smaller changes to (^{ALz)'^^^'^ being 
due to a combination of the weakened spiral amplitude near 
corotation and the change in the value of m. Furthermore, 
we found evidence that the saturation amplitude scales in- 
versely as m, in line with the theory developed by I Sellwood| 
& Binney (2002). Note the simple scaling holds true only 



when the principal dynamical properties, such as Q, X, 
thickness, and gravity softening, are held fixed relative to 
the scale of the mode. Nevertheless, it seems reasonable to 
expect smaller changes to ^(AL^)^^^^^ in general for spirals 

higher m. The exponential decrease of ^(AL^)^^^^^ with in- 
creasing disc thickness applies only for different populations 
subject to the same spiral perturbation. 

We have also run slightly more realistic simulations to 
follow the extent of churning in both thick and thin discs 
that are subject to large numbers of transient spiral waves 
having a variety of rotational symmetries. Fig. [l6] shows 
that changes in the home radii of thick disc particles are 
smaller on average than those of the thin disc, but again the 
tails of the distributions in both populations are almost co- 
extensive. As found in previous work ([Friedli, Benz KeiT 



nicutt|1994||Raboud et a/.|1998||Grenon|1999||Debattista et 



a/.|2006||Minchev et a/.|2011||Bird, Kazantzidis Weinberg 
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2011), the formation of a bar also causes substantial angu- 



lar momentum changes within a disc, but we find that the 
churning effect from multiple spiral patterns still dominates 
changes in the outer disc after the bar has formed. 

We find that vertical action is conserved during radial 
migration, despite the fact that relative changes in our es- 
timated Jz, which can be measured for only about half the 
particles, are as large as ^ 15%. Because the vertical restor- 
ing force to the mid-plane decreases outwards, an increased 
(decreased) home radius causes the particle to experience 
a systematically weaker (stronger) vertical restoring force, 
making it impossible for both vertical energy and action to 
be conserved. We find a clear systematic variation of the 
vertical energy with AL^, but none with AJz leading us 
to conclude that vertical action is the conserved quantity. 
The residual scatter in our measured values of AJz could be 
caused by trapping and escape from multiply-periodic res- 
onant islands in phase space, as well as numerical jitter in 
the iV-body potential, as we discuss at the end of §6. In the 
absence of these complications, we believe that the vertical 
action is conserved. 

Thus conservation of vertical action, and not of vertical 
energy, should be used to prescribe the changes to vertical 
motion in models of chemo-dynamic evolution. It would be 
especially useful to obtain diffusion coefficients that include 
the effects of radial migration as functions of radius and 
height due to a transient spiral or a combination of consecu- 
tive spirals with various corotation radii. We leave this work 
to a future paper. 
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APPENDIX A: NUMERICAL ESTIMATE OF 
VERTICAL ACTION 



Figure A2. Orbits in the meridional plane of the three test par- 
ticles drawn in bold red in the second, fourth and fifth panels 
of Fig. |A1| Note that we show only a small fraction of the full 
orbit used to produce the invariant curves in the (z, z) plane. The 
vertical red lines mark at which the SoS is constructed. The 
initial vertical position and radial and vertical velocities of each 
of these three particles are noted. 
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Fig. |A1| illustrates the SoS for five different energies, all for 
Lz — 7. Each invariant curve is generated by a test particle 
in the initial potential of M2. All the particles in one panel 
have the same energy and angular momentum (noted in the 
top right corner), but differ in the extent of their vertical 
motion. The zero-velocity curve, where a particle's radial 
velocity must be zero for the Lz and E values adopted is 
the outer most green curve. We see that most of phase space 
is regular, but not all consequents mark out simple closed 
invariant curves; some significant fraction of phase space is 
affected by islands in the SoS that surround multiply pe- 
riodic orbits. Chaotic motion, in which consequents fill an 
area rather than lie on a closed curve in the SoS, is extensive 
only for the greatest energy. 

The paths in the meridional plane of three multiply pe- 
riodic orbits are shown in Fig. |A2| They correspond to the 
bold red invariant curves in the second, fourth, and fifth 
SoS panels of Fig. |A1| The last orbit has almost equal ra- 
dial and vertical periods and circulates only clockwise in 
the meridional plane. Since we plot only when the particle 
passes Rh with R > OVb, the consequents lie only in the 
first quadrant of the (z, i) plane. Were we to plot the other 
crossing instead, where R < OVb, the invariant curve would 
be a reflected image in the second quadrant. Also, had we 
chosen the negative root for initial radial velocity, the parti- 
cle would have circulated only counterclockwise and yielded 
a 180°-rotationally symmetric invariant curve in the third 
quadrant. This accounts for the skewness of the invariant 
curves in the SoS, which diminishes for particles of lower 
energy. 

The area in the SoS enclosed by only those orbits having 
a simple invariant curve yields a numerical estimate of the 
vertical action J^QWe find all the consequents for a particle 
as its motion advances over 7 x 10^ dynamical times from 
its position in the frozen, azimuthally averaged, potential. 



and make a numerical estimate of the area integral (eq. 14). 

We repeat this exercise at both the initial and final 
times in the simulation. 



^ A regular orbit trapped about a resonant island has a different 
set of actions that are not of interest here. 
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Figure Al. Surfaces of section in the {z, z) plane, when R = and R > OVq, for test particles in the initial potential of simulation M2. 
All particles in each panel have the same total energy and corresponding Zrmmax, which are noted in the top right corner. The outer 
most green curves mark the zero-velocity curves, at which R = OVq for the given z, Lz and E. Finally, the bold red curves are from the 
three orbits shown in Fig. |A2| 
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